Density matrix renormalization group (DMRG) is one of the most important numerical methods for low dimensional strongly correlated lattice models.^[1–3] The ground state of a one-dimensional (1D) system can be obtained variationally by DMRG with extremely high efficiency and precision. After its invention for the real space lattice model, the application of DMRG has been extended to the momentum space^[4,5] and quantum chemistry.^[6–8] Besides the ground state algorithm, a serial of techniques, including the time-dependent and finite-temperature DMRG, are generalized and developed.^[9–15] In addition, many novel quantum phases in two-dimensional (2D) interaction systems have been investigated by DMRG, e.g., stripe phase in high temperature superconductor^[16–22] and quantum spin liquid.^[23–25] However, in these extended applications, especially for the two and quasi-two dimensional problems, the computational cost of DMRG increases dramatically for the convergent results. These lead to great challenges to the computing and storage capabilities of the contemporary computers.

Nowadays, variety of optimized methods have been proposed to reduce the computational cost of DMRG, for example, targeting the states in smaller Hilbert subspace in the presence of multiple symmetries,^[5,26,27] the dynamical block state selection approach,^[28,29] one-site sweep strategy,^[30,31] and good initial wave function for the eigensolver.^[32] On the other hand, to take full advantages of the advanced high-performance computing technique and further shorten the time cost of DMRG, there have been different parallelization strategies on high-performance computers, e.g., real space parallel DMRG,^[33] shared/distributed-memory multi-core parallelization,^[34,35] and parallelization on heterogeneous architectures.^[36,37] Especially, the heterogeneous architectures employing both central processing unit (CPU) and graphics processing unit (GPU) have potential to significantly improve the efficiency, since the floating-point performance and memory bandwidth of GPU are usually many times higher than the host part.

Along this direction, a pioneer work has implemented DMRG on the hybrid CPU–GPU architecture, and achieved considerable speedup.^[36] However, this implementation is based on the two-subblock configuration, the memory cost of which is d² (d is the local Hilbert space dimension of a single site) times larger than that of the four-subblock expression. Recently, we propose a new hybrid CPU–GPU parallel strategy with the Hamiltonian expressed by the operators of four subblocks, which greatly reduces the GPU memory cost.^[37] In our implementation, one is able to use much larger number of DMRG kept states (upto 10⁴) to investigate more complex systems, such as spinfull fermions on the quasi-two-dimensional lattice.

In our previous work, the hybrid strategy focuses on the diagonalization of the Hamiltonian, which is the most time-consuming part for the DMRG procedure with large number of kept states. According to a dynamical load balancing strategy, the matrix and vector operations in diagonalizing the Hamiltonian are distributed to CPU and GPU, and are performed using parallel MKL^[38] and CUBLAS^[39] subroutines, respectively. While these libraries generally achieve high performance in most cases, they are not always the best choice. For example, we get low performance when using these parallel subroutines for small matrix and vector operations. Therefore, in this work we improve our parallelization strategy by optimizing the diagonalization part. On the host, when acting the Hamiltonian on the wavefunction, we adopt OpenMP application programming interface^[40] and sequential MKL subroutines, instead of parallel MKL subroutines. On the GPU, for vector operations in acting the Hamiltonian on the wavefunction, we use our own subroutines which gets higher bandwidth for linear combination of vectors. In addition, the previous DMRG truncation procedure with lots of matrix operations is only performed by the CPU. In the current work, we perform most of these operations on the GPU to further reduce the total computational time.

This paper is organized as follows. In Section 2, we review the finite-system DMRG and our previous parallel strategy, and get the time cost for each part in calculating the ground state of the 4-leg Hubbard model with next-nearest hopping. In Section 3, we introduce specifically the optimized hybrid parallel strategy for two computing sections, the diagonalization of the Hamiltonian and the truncation procedure. In Section 4, we show the numerical results and the performance benchmark comparing with the previous work. The conclusion is given in Section 5.

A DMRG routine for finite lattice usually contains 1) infinite-system DMRG with a small number of states to produce a good initial wavefunction, and 2) a few finite-system DMRG sweeps to get the convergent result. For the quasi-two-dimensional lattice models, finite DMRG sweeps take most of the computational time. In each sweep, the wavefunction is updated through the lattice system from left to right, and then back to the left side. Without loss of generality, we briefly introduce one DMRG step with sweeping through the lattice system from left to right as an example. In our implementation, the lattice system is partitioned into four subblocks (see Fig. 1): block S, block E, site m₁, and site m₂.

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. Schematic one-step procedure of finite-system DMRG sweeps in four-subblock formulation.

A wave function |ψ〉 can be expressed in four-subblock basis as

To demonstrate the performance of our improved implementation, we apply it to get the ground state of the 4-leg Fermi Hubbard ladder with the next-nearest hopping, which is defined as

In our implementation, we employ these two good quantum numbers to reduces the computational and memory cost. In this paper, we set the nearest hopping t = 1 as the energy scale, and fix the next-nearest hopping t′ = 0.25 and the Coulomb repulsion U = 6. For the performance benchmark, we focus on the 4-leg ladder with leg-length L_x = 16 and particle filling 7/8. The cylinder geometry is adopted, with open and periodic boundaries along the leg and rung, respectively. For all calculations, we take 5 sweeps to obtain the convergent energy.

Employing the previous hybrid parallel strategy^[37] on the next-nearest hopping Hubbard ladder, we get the computational time proportion of several computing sections, as shown in Fig. 2. In the previous parallel scheme, applying the Hamiltonian to the wavefunction takes more than 60% of the total time cost. On the CPU, totally relying on MKL subroutines may not take good advantage of the computational resources for the small scale matrix and vector operations. On the GPU, we get low memory bandwidth for linear combination of vectors with CUBLAS subroutines based on Computer-Unified Device Architecture (CUDA). The time cost ratio of the truncation procedure is about 20%, which contains diagonalizing the reduced density matrix and projecting the operator on the new basis. Applying the Hamiltonian to the wavefunction and the truncation procedure takes most of the computational time (about 80% to 85%), and these two parts still have space to be improved. In the next section, we introduce explicitly how we further optimize the hybrid parallelization of the DMRG program with a large number of kept states.

Fig. 2.

Figure Option ViewDownloadNew Window

Fig. 2. The computational time ratio to the total time cost for the most time-consuming parts. Applying the Hamiltonian to the wavefunction (|ϕ〉 = H|ψ〉, denoted by ×) is the core operation for diagonalizing the Hamiltonian. The truncation procedure (°) contains diagonalizing the reduced density matrix (▽) and projecting the operator on the new basis (Δ). The former can be considered as the singular value decomposition (SVD) of the superblock wavefunction.

To test the performance of our optimized parallel strategy, we apply it to calculate the ground state of the model (3), with different DMRG kept states (D ∈ {4096, 6144, 8192, 10240}). In our hybrid parallel environment, there are 1 CPU (Intel Xeon E5-2650v3, 10 cores, 2.3 GHz) and 1 GPU (12 GB GPU memory) of the NVIDIA Tesla K80 graphic card. The maximum number of GPU kernels executed concurrently in our calculations is 16. To make sure our numerical code is correct, we compare our results with the DMRG method in Itensor^[43] by calculating the energy per site of the ground state. As shown in Table 1, these two implementations give the accordance results.

Table 1.

Table 1.

Table 1. The ground state energy per site E0/L. L is the total site number of the lattice. CPU–GPU denotes our results. . D CPU–GPU Itensor 4096 −0.9280790 −0.9280679 6144 −0.9284130 −0.9284069 8192 −0.9285577 −0.9285539 10240 −0.9286339 −0.9286315

Table 1. The ground state energy per site E0/L. L is the total site number of the lattice. CPU–GPU denotes our results. .

We also compute the charge density of the ground state, and display the hole density in Fig. 4. As addressed in Ref. [22], the ground state of the model with fixed parameters (see Section 2) displays a stripe phase with charge density waves, which is usually connected to high temperature superconductors. As shown in Fig. 4, the locations of the maximum hole densities, which are also believed to be the locations of the magnetic domain walls, divide the real space into unit-cells with wavelength λ = 4. These results are in agreement with studies in Ref. [22].

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. The hole density on each site. Hole density is represented by both the area and the color depth of the disk. Black solid lines denote hopping, and red dashed lines represent the positions of maximum hole density. Here D = 10240.

For different DMRG kept states D, we set the total computational time of our previous strategy as the unit, and demonstrate the performance improvement of each optimized section by showing the relative computational time reduction, as displayed in Fig. 5. In acting the Hamiltonian on the wavefunction, we can save more than 20% of the total time cost for D = 4096 by optimizing operations on both CPU and GPU. The performance of parallel MKL subroutines we used before on the CPU has better performance for larger scale operations, so that the improvement in this part looks more prominent for smaller kept states. In the truncation procedure, the matrix operations always achieve high performance on GPU than CPU for large scale matrices. Therefore, we further take advantage of the GPU in our new parallel scheme to save more computational time. In this part, the performance improvement increases as D increases. The total computational time of our improved parallel strategy is about 23%–29% less than that in the previous work.

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. The relative time cost reduction of different computing sections compared to the total computational time in our previous work. Truncation procedure (°) contains SVD (▽) and projection (Δ).

In our hybrid parallel strategy, we have considerably improved the performance of matrix and vector operations in DMRG calculation of the quasi-2D Hubbard model. Meanwhile, in the procedure of acting the Hamiltonian on the wave function (see Algorithm 1), which is the most time-consuming part, one may achieve better performance in the following aspects. First, for some models, there might be many small subspaces involved, which can lead to low performance in steps 1 and 3. In these cases, we can apply high performance subroutines specially for small matrix operations. Second, high host memory bandwidth may be achieved by extending the strategy optimizing the operations of step 2 in GPU to the host part as well. Third, if the floating-point operations fall in a small number of subspaces, it is always difficult to balance the load between CPU and GPU. At this point, we can divide the states in large subspaces into many parts, and distribute the according operations into many smaller tasks.

In this paper, based on our previous work, we propose an optimized parallel scheme of DMRG method on the hybrid CPU–GPU architecture. Applying it to solve the ground state of a 4-leg Hubbard model with next-nearest hopping, we obtain the charge stripes usually being observed in high temperature superconductors, and demonstrate the performance improvement. Compared to the previous strategy, both the diagonalization and the truncation procedure are optimized, and these optimizations significantly reduce the total computational time in the DMRG calculation. Our results indicate the advantage of the new hybrid parallel strategy and the importance of high performance computing in strongly correlated algorithms. We expect our hybrid parallel strategy benefits DMRG applications with large number of kept states, for example, ground state DMRG for quasi-2D lattice models and time-dependent DMRG methods.